knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = TRUE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE # ! set to TRUE !

if (get_bayesfactor) {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = paste0('_final_with_plan_', as.character(iterations))
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1635 0.0000 5.0000 241

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Own Actionplan",
    'Partner Actionplan',
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,11),
  "Between-Person Effects" = c(12,18),
  "Random Effects" = c(19, 25), 
  "Additional Parameters" = c(26,26)
  )


rows_to_pack_ordinal <- list(
  "Within-Person Effects" = c(2+5,11+5),
  "Between-Person Effects" = c(12+5,18+5),
  "Random Effects" = c(19+5, 25+5), 
  "Additional Parameters" = c(26+5,26+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
if (check_models) {
  check_brms(pa_sub, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')
}
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 7, observations = 3736, p-value = 0.56
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000267666 0.002817184
## sample estimates:
## outlier frequency (expected: 0.00149892933618844 ) 
##                                        0.001873662
if (do_priorsense) {
  priorsense_vars <- c(
      'Intercept',
      'b_persuasion_self_cw',
      'b_persuasion_partner_cw',
      'b_pressure_self_cw',
      'b_pressure_partner_cw',
      'b_pushing_self_cw',
      'b_pushing_partner_cw'
  )
  
  hurdle_priorsense_vars <- c(
    'Intercept_hu',
    'b_hu_persuasion_self_cw',
    'b_hu_persuasion_partner_cw',
    'b_hu_pressure_self_cw',
    'b_hu_pressure_partner_cw',
    'b_hu_pushing_self_cw',
    'b_hu_pushing_partner_cw'
  )
  
  gc()
  priorsense::powerscale_sensitivity(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_dens(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_ecdf(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_quantities(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
}
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, stats_to_report = stats_to_report, rope_range
## = rope_range_continuous, : Coefficients were exponentiated. Double check if
## this was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
Intercept 0.26*** 0.04 [ 0.19, 0.35] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 10801 19932 37.92*** 2.64 [33.02, 43.54] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 12653 22304
Within-Person Effects
Daily persuasion experienced 1.59*** 0.11 [ 1.39, 1.85] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 26190 27004 1.03 0.03 [ 0.98, 1.09] 0.859 [0.92, 1.08] 0.967 0.019 Very Strong Evidence for Null 1.000 17362 24808
Daily persuasion utilized (partner’s view) 1.34*** 0.09 [ 1.19, 1.54] 1.000 [0.84, 1.20] 0.035 >100 Overwhelming Evidence 1.000 30919 27774 1.03 0.02 [ 0.99, 1.08] 0.916 [0.92, 1.08] 0.978 0.023 Very Strong Evidence for Null 1.000 22959 28221
Daily pressure experienced 0.95 0.14 [ 0.69, 1.29] 0.622 [0.84, 1.20] 0.742 0.079 Strong Evidence for Null 1.000 35225 27790 0.91* 0.04 [ 0.82, 1.00] 0.975 [0.92, 1.08] 0.343 0.154 Moderate Evidence for Null 1.000 36639 29172
Daily pressure utilized (partner’s view) 1.49* 0.28 [ 1.05, 2.29] 0.987 [0.84, 1.20] 0.113 1.098 Weak Evidence 1.000 34454 22964 0.95 0.04 [ 0.86, 1.03] 0.894 [0.92, 1.08] 0.700 0.038 Strong Evidence for Null 1.000 37501 28659
Daily pushing experienced 0.94 0.14 [ 0.71, 1.28] 0.665 [0.84, 1.20] 0.735 0.083 Strong Evidence for Null 1.000 23390 25680 0.99 0.03 [ 0.92, 1.05] 0.674 [0.92, 1.08] 0.966 0.014 Very Strong Evidence for Null 1.000 31510 30242
Daily pushing utilized (partner’s view) 1.28* 0.15 [ 1.04, 1.64] 0.990 [0.84, 1.20] 0.264 0.819 Weak Evidence for Null 1.000 35469 27582 0.96 0.03 [ 0.90, 1.02] 0.914 [0.92, 1.08] 0.882 0.031 Strong Evidence for Null 1.000 33140 30422
Day 0.86 0.13 [ 0.64, 1.15] 0.843 [0.84, 1.20] 0.566 0.124 Moderate Evidence for Null 1.000 61855 30632 0.99 0.06 [ 0.88, 1.11] 0.577 [0.92, 1.08] 0.790 0.025 Very Strong Evidence for Null 1.000 50335 29874
Own Actionplan 9.48*** 0.97 [ 7.77, 11.60] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 45197 31239 1.33*** 0.06 [ 1.21, 1.45] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 46130 29477
Partner Actionplan 1.17 0.12 [ 0.96, 1.42] 0.941 [0.84, 1.20] 0.583 0.178 Moderate Evidence for Null 1.000 44112 30129 1.08 0.05 [ 0.99, 1.17] 0.962 [0.92, 1.08] 0.518 0.084 Strong Evidence for Null 1.000 47984 32207
Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.40 0.50 [ 0.69, 2.84] 0.825 [0.84, 1.20] 0.256 0.280 Moderate Evidence for Null 1.001 8405 15769 1.03 0.16 [ 0.77, 1.40] 0.590 [0.92, 1.08] 0.386 0.061 Strong Evidence for Null 1.000 7638 14348
Mean persuasion utilized (partner’s view) 1.33 0.47 [ 0.65, 2.67] 0.783 [0.84, 1.20] 0.289 0.240 Moderate Evidence for Null 1.000 8466 14879 0.99 0.15 [ 0.73, 1.34] 0.524 [0.92, 1.08] 0.396 0.061 Strong Evidence for Null 1.000 7814 14204
Mean pressure experienced 0.40* 0.17 [ 0.18, 0.90] 0.986 [0.84, 1.20] 0.034 2.480 Weak Evidence 1.000 14450 24340 1.18 0.21 [ 0.83, 1.67] 0.826 [0.92, 1.08] 0.227 0.110 Moderate Evidence for Null 1.000 11628 20761
Mean pressure utilized (partner’s view) 0.51 0.21 [ 0.22, 1.15] 0.948 [0.84, 1.20] 0.096 0.776 Weak Evidence for Null 1.000 14670 22245 0.92 0.17 [ 0.65, 1.32] 0.671 [0.92, 1.08] 0.300 0.081 Strong Evidence for Null 1.000 12201 21614
Mean pushing experienced 1.17 0.62 [ 0.41, 3.28] 0.616 [0.84, 1.20] 0.252 0.278 Moderate Evidence for Null 1.000 12680 20855 1.21 0.27 [ 0.77, 1.88] 0.797 [0.92, 1.08] 0.193 0.128 Moderate Evidence for Null 1.000 11020 18100
Mean pushing utilized (partner’s view) 2.07 1.10 [ 0.71, 5.84] 0.916 [0.84, 1.20] 0.105 0.687 Weak Evidence for Null 1.000 12400 20808 1.30 0.30 [ 0.82, 2.06] 0.871 [0.92, 1.08] 0.141 0.175 Moderate Evidence for Null 1.000 11072 19999
Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.73 0.10 [0.56, 0.99] NA NA NA NA NA 1.000 10861 19564 0.30 0.04 [0.23, 0.41] NA NA NA NA NA 1.000 10982 19107
sd(Daily persuasion experienced) 0.21 0.09 [0.03, 0.40] NA NA NA NA NA 1.001 10174 8865 0.11 0.02 [0.07, 0.17] NA NA NA NA NA 1.000 20621 28034
sd(Daily persuasion utilized (partner’s view)) 0.19 0.09 [0.02, 0.37] NA NA NA NA NA 1.000 9864 8686 0.08 0.02 [0.04, 0.13] NA NA NA NA NA 1.000 21526 23726
sd(Daily pressure experienced) 0.18 0.16 [0.01, 0.68] NA NA NA NA NA 1.000 16735 17846 0.06 0.05 [0.00, 0.22] NA NA NA NA NA 1.000 16670 18381
sd(Daily pressure utilized (partner’s view)) 0.23 0.21 [0.01, 0.89] NA NA NA NA NA 1.000 16584 19362 0.06 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 17699 15819
sd(Daily pushing experienced) 0.57 0.16 [0.30, 0.98] NA NA NA NA NA 1.000 18443 25123 0.09 0.04 [0.01, 0.17] NA NA NA NA NA 1.000 11617 9522
sd(Daily pushing utilized (partner’s view)) 0.22 0.14 [0.02, 0.54] NA NA NA NA NA 1.000 13769 12571 0.08 0.03 [0.01, 0.15] NA NA NA NA NA 1.000 12829 10346
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 0.68 0.01 [0.65, 0.70] NA NA NA NA NA 1.000 54627 28168 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.76), b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.74). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

Hurdle part of the model on the left, non-zero part towards the right side of the table

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## This is posterior version 1.6.0
## 
## Attaching package: 'posterior'
## The following object is masked from 'package:bayesplot':
## 
##     rhat
## The following objects are masked from 'package:stats':
## 
##     mad, sd, var
## The following objects are masked from 'package:base':
## 
##     %in%, match
## Registering fonts with R
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 80 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00828

$persuasion_partner_cw

## Warning: Removed 61 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00653

$pressure_self_cw

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.011

$pressure_partner_cw

## Picking joint bandwidth of 0.0184

$pushing_self_cw

## Picking joint bandwidth of 0.0103

$pushing_partner_cw

## Picking joint bandwidth of 0.0101

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  )

home_dir <- getwd()
output_dir <- file.path(home_dir, 'Output', 'Plots')

for (i in 1:length(conds_eff)) {
  effname <- names(conds_eff)[i]
  eff_plot <- conds_eff[[i]]
  x_label_i <- x_label[[i]]
  
  rmarkdown::render(
    file.path(output_dir, 'BeautifulPlotWithNote.Rmd'), 
    output_file = file.path(output_dir, paste0('Graphic_', effname, '.pdf')),
    params = list(
      home_dir = home_dir,
      output_dir = output_dir,
      p_i = eff_plot,
      p_name = effname,
      x_label = x_label_i
      ),
    envir = new.env(),
    quiet = TRUE
  )
}

print('done')

Comparing effect size of pressure and pushing

hypothesis(pa_sub, "pressure_self_cw < pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... < 0    -0.08      0.06    -0.19     0.02       9.59
##   Post.Prob Star
## 1      0.91     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
if (check_models) {
  check_brms(pa_obj_log, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
}
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 28, observations = 3337, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.001198681 0.004794726
## sample estimates:
## outlier frequency (expected: 0.00297872340425532 ) 
##                                         0.00839077
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(pa_obj_log, variable = priorsense_vars)
}
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, stats_to_report = stats_to_report, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 111.17*** 5.93 [99.95, 123.93] 1.000 [0.94, 1.07] 0.000 >100 Overwhelming Evidence 1.001 6903 13661
Within-Person Effects
Daily persuasion experienced 1.03 0.02 [ 1.00, 1.06] 0.953 [0.94, 1.07] 0.995 0.026 Very Strong Evidence for Null 1.000 26993 29364
Daily persuasion utilized (partner’s view) 1.02 0.02 [ 0.99, 1.05] 0.864 [0.94, 1.07] 0.998 0.011 Very Strong Evidence for Null 1.000 31647 30752
Daily pressure experienced 0.95 0.03 [ 0.88, 1.01] 0.947 [0.94, 1.07] 0.602 0.055 Strong Evidence for Null 1.000 47423 30997
Daily pressure utilized (partner’s view) 0.98 0.03 [ 0.92, 1.05] 0.694 [0.94, 1.07] 0.916 0.015 Very Strong Evidence for Null 1.000 50157 31658
Daily pushing experienced 1.02 0.02 [ 0.96, 1.07] 0.731 [0.94, 1.07] 0.974 0.012 Very Strong Evidence for Null 1.000 35878 29732
Daily pushing utilized (partner’s view) 1.00 0.02 [ 0.96, 1.04] 0.510 [0.94, 1.07] 0.997 0.008 Very Strong Evidence for Null 1.000 48187 32131
Day 0.97 0.03 [ 0.91, 1.04] 0.788 [0.94, 1.07] 0.845 0.020 Very Strong Evidence for Null 1.000 77038 29175
Own Actionplan 1.06* 0.03 [ 1.01, 1.12] 0.993 [0.94, 1.07] 0.530 0.225 Moderate Evidence for Null 1.000 57658 31464
Partner Actionplan 1.05 0.03 [ 1.00, 1.10] 0.972 [0.94, 1.07] 0.753 0.060 Strong Evidence for Null 1.000 57893 31743
Daily weartime 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.07] 1.000 0.053 Strong Evidence for Null 1.000 41154 26326
Between-Person Effects
Mean persuasion experienced 1.10 0.16 [ 0.83, 1.46] 0.756 [0.94, 1.07] 0.279 0.072 Strong Evidence for Null 1.001 5287 11364
Mean persuasion utilized (partner’s view) 0.98 0.14 [ 0.73, 1.30] 0.562 [0.94, 1.07] 0.344 0.059 Strong Evidence for Null 1.001 5210 10807
Mean pressure experienced 0.99 0.14 [ 0.74, 1.33] 0.531 [0.94, 1.07] 0.344 0.058 Strong Evidence for Null 1.001 7594 15835
Mean pressure utilized (partner’s view) 0.97 0.14 [ 0.74, 1.29] 0.577 [0.94, 1.07] 0.347 0.058 Strong Evidence for Null 1.001 6847 13727
Mean pushing experienced 0.95 0.19 [ 0.63, 1.43] 0.600 [0.94, 1.07] 0.241 0.084 Strong Evidence for Null 1.001 8398 15531
Mean pushing utilized (partner’s view) 1.22 0.25 [ 0.82, 1.84] 0.840 [0.94, 1.07] 0.156 0.133 Moderate Evidence for Null 1.001 8149 15791
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.916 [0.94, 1.07] 1.000 0.000 Very Strong Evidence for Null 1.000 50144 33445
Random Effects
sd(Intercept) 0.29 0.04 [0.23, 0.39] NA NA NA NA NA 1.001 8384 17687
sd(Daily persuasion experienced) 0.05 0.01 [0.02, 0.08] NA NA NA NA NA 1.000 16704 10988
sd(Daily persuasion utilized (partner’s view)) 0.05 0.02 [0.03, 0.09] NA NA NA NA NA 1.000 19719 20938
sd(Daily pressure experienced) 0.04 0.03 [0.00, 0.13] NA NA NA NA NA 1.000 19449 21382
sd(Daily pressure utilized (partner’s view)) 0.03 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 23291 19151
sd(Daily pushing experienced) 0.07 0.04 [0.01, 0.15] NA NA NA NA NA 1.000 9926 12203
sd(Daily pushing utilized (partner’s view)) 0.03 0.03 [0.00, 0.09] NA NA NA NA NA 1.000 16153 18099
Additional Parameters
sigma 0.57 0.01 [0.56, 0.59] NA NA NA NA NA 1.000 65401 28026 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.89), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.74), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.74), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.72), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.77), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.82). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
if (check_models) {
  check_brms(mood_gauss, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss, integer = FALSE)
}
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 30, observations = 3736, p-value =
## 4.114e-10
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.005424162 0.011443600
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.008029979
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(mood_gauss, variable = priorsense_vars)
}
summarize_brms(
  mood_gauss, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.67*** 0.10 [ 3.46, 3.88] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.001 4513 9768
Within-Person Effects
Daily persuasion experienced 0.00 0.02 [-0.04, 0.04] 0.523 [-0.11, 0.11] 1.000 0.004 Very Strong Evidence for Null 1.000 35272 26862
Daily persuasion utilized (partner’s view) 0.02 0.02 [-0.02, 0.07] 0.810 [-0.11, 0.11] 1.000 0.007 Very Strong Evidence for Null 1.000 30965 30203
Daily pressure experienced -0.03 0.05 [-0.14, 0.07] 0.717 [-0.11, 0.11] 0.937 0.012 Very Strong Evidence for Null 1.000 36650 27776
Daily pressure utilized (partner’s view) -0.03 0.05 [-0.15, 0.08] 0.694 [-0.11, 0.11] 0.930 0.012 Very Strong Evidence for Null 1.000 34035 26217
Daily pushing experienced 0.01 0.03 [-0.06, 0.07] 0.577 [-0.11, 0.11] 0.999 0.007 Very Strong Evidence for Null 1.000 42746 31073
Daily pushing utilized (partner’s view) 0.07* 0.03 [ 0.00, 0.14] 0.976 [-0.11, 0.11] 0.895 0.058 Strong Evidence for Null 1.000 31862 28826
Day 0.26*** 0.06 [ 0.15, 0.37] 1.000 [-0.11, 0.11] 0.005 >100 Overwhelming Evidence 1.000 58650 29078
Own Actionplan 0.10** 0.04 [ 0.03, 0.18] 0.996 [-0.11, 0.11] 0.618 0.268 Moderate Evidence for Null 1.000 48149 31839
Partner Actionplan -0.03 0.04 [-0.11, 0.05] 0.792 [-0.11, 0.11] 0.982 0.011 Very Strong Evidence for Null 1.000 47941 31849
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.34 0.27 [-0.20, 0.90] 0.892 [-0.11, 0.11] 0.154 0.120 Moderate Evidence for Null 1.001 3933 7303
Mean persuasion utilized (partner’s view) 0.22 0.27 [-0.32, 0.78] 0.794 [-0.11, 0.11] 0.236 0.076 Strong Evidence for Null 1.001 3906 7095
Mean pressure experienced -0.30 0.27 [-0.85, 0.24] 0.866 [-0.11, 0.11] 0.181 0.102 Moderate Evidence for Null 1.001 4698 9345
Mean pressure utilized (partner’s view) -0.31 0.27 [-0.86, 0.22] 0.877 [-0.11, 0.11] 0.170 0.110 Moderate Evidence for Null 1.001 4686 9233
Mean pushing experienced 0.21 0.38 [-0.55, 0.97] 0.707 [-0.11, 0.11] 0.204 0.087 Strong Evidence for Null 1.000 6174 12395
Mean pushing utilized (partner’s view) 0.37 0.38 [-0.39, 1.13] 0.835 [-0.11, 0.11] 0.149 0.124 Moderate Evidence for Null 1.000 6234 12509
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.60 0.07 [0.48, 0.78] NA NA NA NA NA 1.000 7398 14380
sd(Daily persuasion experienced) 0.04 0.03 [0.00, 0.10] NA NA NA NA NA 1.000 10459 14097
sd(Daily persuasion utilized (partner’s view)) 0.07 0.03 [0.01, 0.13] NA NA NA NA NA 1.001 11847 10659
sd(Daily pressure experienced) 0.07 0.06 [0.00, 0.24] NA NA NA NA NA 1.000 16312 19973
sd(Daily pressure utilized (partner’s view)) 0.08 0.07 [0.00, 0.26] NA NA NA NA NA 1.000 15919 18986
sd(Daily pushing experienced) 0.05 0.04 [0.00, 0.15] NA NA NA NA NA 1.000 15239 17528
sd(Daily pushing utilized (partner’s view)) 0.07 0.05 [0.00, 0.17] NA NA NA NA NA 1.000 13634 14941
Additional Parameters
sigma 0.96 0.01 [0.94, 0.98] NA NA NA NA NA 1.000 54740 29154 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.82), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.8), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.8), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.83), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.77), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.89). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
if (check_models) {
  check_brms(reactance_ordinal)
  DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')
}
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 4 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 2, observations = 756, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.001322751
## sample estimates:
## outlier frequency (expected: 0.000317460317460317 ) 
##                                         0.002645503
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(reactance_ordinal, variable = priorsense_vars)
}
summarize_brms(
  reactance_ordinal, 
  stats_to_report = stats_to_report,
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
## Sampling priors, please wait...
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 3.36*** 1.02 [ 1.88, 6.23] 1.000 [0.84, 1.20] 0.000 69.572 Very Strong Evidence 1.000 33161 30736
Intercept[2] 7.31*** 2.30 [ 4.00, 13.82] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 33914 29999
Intercept[3] 20.38*** 6.87 [ 10.74, 40.28] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 35591 31079
Intercept[4] 89.21*** 34.75 [ 42.90, 195.22] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 39807 32159
Intercept[5] 3029.60*** 1996.62 [916.76, 12308.87] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 49408 30767
Within-Person Effects
Daily persuasion experienced 0.84* 0.07 [ 0.71, 0.99] 0.982 [0.84, 1.20] 0.549 0.297 Moderate Evidence for Null 1.000 33919 28293
Daily persuasion utilized (partner’s view) 1.02 0.10 [ 0.83, 1.23] 0.583 [0.84, 1.20] 0.925 0.040 Strong Evidence for Null 1.000 32260 28470
Daily pressure experienced 1.84* 0.37 [ 1.16, 2.66] 0.992 [0.84, 1.20] 0.031 2.539 Weak Evidence 1.000 21582 22822
Daily pressure utilized (partner’s view) 1.24 0.30 [ 0.69, 2.11] 0.807 [0.84, 1.20] 0.376 0.158 Moderate Evidence for Null 1.000 23305 19657
Daily pushing experienced 1.22 0.13 [ 0.98, 1.52] 0.962 [0.84, 1.20] 0.445 0.217 Moderate Evidence for Null 1.000 28606 27455
Daily pushing utilized (partner’s view) 0.94 0.12 [ 0.71, 1.22] 0.693 [0.84, 1.20] 0.776 0.059 Strong Evidence for Null 1.000 32445 27556
Day 1.48 0.51 [ 0.75, 2.88] 0.872 [0.84, 1.20] 0.221 0.266 Moderate Evidence for Null 1.000 47791 31765
Own Actionplan 0.85 0.24 [ 0.49, 1.49] 0.714 [0.84, 1.20] 0.410 0.134 Moderate Evidence for Null 1.000 35305 32306
Partner Actionplan 0.92 0.24 [ 0.55, 1.52] 0.632 [0.84, 1.20] 0.493 0.110 Moderate Evidence for Null 1.000 37423 31789
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.10 0.56 [ 0.39, 3.02] 0.572 [0.84, 1.20] 0.269 0.211 Moderate Evidence for Null 1.000 15958 22898
Mean persuasion utilized (partner’s view) 1.37 0.77 [ 0.45, 4.27] 0.709 [0.84, 1.20] 0.219 0.265 Moderate Evidence for Null 1.000 16443 22952
Mean pressure experienced 3.52* 1.93 [ 1.22, 10.68] 0.990 [0.84, 1.20] 0.020 3.211 Moderate Evidence 1.000 18630 25099
Mean pressure utilized (partner’s view) 1.16 0.67 [ 0.36, 3.54] 0.602 [0.84, 1.20] 0.235 0.236 Moderate Evidence for Null 1.000 17309 24737
Mean pushing experienced 1.25 0.93 [ 0.30, 5.72] 0.623 [0.84, 1.20] 0.187 0.310 Weak Evidence for Null 1.000 20161 26058
Mean pushing utilized (partner’s view) 0.11* 0.10 [ 0.02, 0.65] 0.993 [0.84, 1.20] 0.009 7.349 Moderate Evidence 1.000 26207 29097
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.80 0.20 [0.46, 1.27] NA NA NA NA NA 1.000 15427 22186
sd(Daily persuasion experienced) 0.16 0.12 [0.01, 0.43] NA NA NA NA NA 1.000 7784 15596
sd(Daily persuasion utilized (partner’s view)) 0.21 0.14 [0.01, 0.52] NA NA NA NA NA 1.001 10675 13886
sd(Daily pressure experienced) 0.56 0.25 [0.11, 1.18] NA NA NA NA NA 1.000 8950 7907
sd(Daily pressure utilized (partner’s view)) 0.42 0.40 [0.02, 1.56] NA NA NA NA NA 1.001 9732 18185
sd(Daily pushing experienced) 0.21 0.14 [0.01, 0.51] NA NA NA NA NA 1.000 11006 14481
sd(Daily pushing utilized (partner’s view)) 0.15 0.14 [0.01, 0.62] NA NA NA NA NA 1.000 17161 19837
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.8), b_Intercept[4] and b_Intercept[3] (r = 0.86), b_pressure_self_cb
##   and b_persuasion_self_cb (r = 0.71), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.79). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="01_FinalModelsPlan_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
if (check_models) {
  check_brms(is_reactance)
  DHARMa.check_brms.all(is_reactance, integer = FALSE)
}
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 35 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 2, observations = 756, p-value = 0.6663
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0003205434 0.0095234999
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.002645503
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(is_reactance, variable = priorsense_vars)
}
summarize_brms(
  is_reactance, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.34*** 0.11 [0.17, 0.64] 1.000 [0.83, 1.20] 0.003 12.437 Strong Evidence 1.000 28741 31565
Within-Person Effects
Daily persuasion experienced 0.84 0.08 [0.68, 1.01] 0.968 [0.83, 1.20] 0.527 0.215 Moderate Evidence for Null 1.000 26792 22773
Daily persuasion utilized (partner’s view) 1.12 0.16 [0.84, 1.53] 0.780 [0.83, 1.20] 0.666 0.074 Strong Evidence for Null 1.000 20781 21793
Daily pressure experienced 1.97* 0.64 [1.01, 4.43] 0.976 [0.83, 1.20] 0.053 1.201 Weak Evidence 1.000 16708 16138
Daily pressure utilized (partner’s view) 1.42 0.60 [0.58, 4.18] 0.804 [0.83, 1.20] 0.234 0.241 Moderate Evidence for Null 1.000 16914 17573
Daily pushing experienced 1.33* 0.18 [1.03, 1.75] 0.986 [0.83, 1.20] 0.208 0.603 Weak Evidence for Null 1.000 24947 25401
Daily pushing utilized (partner’s view) 0.92 0.17 [0.62, 1.36] 0.677 [0.83, 1.20] 0.605 0.087 Strong Evidence for Null 1.000 28545 25800
Day 1.66 0.65 [0.77, 3.60] 0.905 [0.83, 1.20] 0.160 0.368 Weak Evidence for Null 1.000 39966 30260
Own Actionplan 0.87 0.28 [0.46, 1.64] 0.665 [0.83, 1.20] 0.395 0.139 Moderate Evidence for Null 1.000 30891 30274
Partner Actionplan 0.87 0.26 [0.49, 1.54] 0.681 [0.83, 1.20] 0.421 0.130 Moderate Evidence for Null 1.000 35248 29799
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.90 1.16 [0.57, 6.60] 0.856 [0.83, 1.20] 0.135 0.430 Weak Evidence for Null 1.000 14375 21061
Mean persuasion utilized (partner’s view) 1.83 1.23 [0.50, 7.20] 0.819 [0.83, 1.20] 0.144 0.395 Weak Evidence for Null 1.000 14643 22514
Mean pressure experienced 18.34** 19.65 [2.51, 165.74] 0.998 [0.83, 1.20] 0.002 29.292 Strong Evidence 1.000 18604 23442
Mean pressure utilized (partner’s view) 2.29 2.54 [0.23, 19.25] 0.769 [0.83, 1.20] 0.095 0.613 Weak Evidence for Null 1.000 17099 24008
Mean pushing experienced 0.86 0.85 [0.12, 6.20] 0.562 [0.83, 1.20] 0.144 0.401 Weak Evidence for Null 1.000 16631 23514
Mean pushing utilized (partner’s view) 0.08* 0.09 [0.01, 0.68] 0.990 [0.83, 1.20] 0.009 6.280 Moderate Evidence 1.000 21332 27410
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.15 0.25 [0.74, 1.76] NA NA NA NA NA 1.000 12310 20831
sd(Daily persuasion experienced) 0.22 0.14 [0.01, 0.51] NA NA NA NA NA 1.000 7209 12307
sd(Daily persuasion utilized (partner’s view)) 0.49 0.20 [0.13, 0.98] NA NA NA NA NA 1.000 11110 10788
sd(Daily pressure experienced) 1.08 0.55 [0.14, 2.43] NA NA NA NA NA 1.001 6642 6944
sd(Daily pressure utilized (partner’s view)) 0.83 0.68 [0.04, 2.75] NA NA NA NA NA 1.000 10859 14383
sd(Daily pushing experienced) 0.25 0.16 [0.02, 0.62] NA NA NA NA NA 1.000 10817 11169
sd(Daily pushing utilized (partner’s view)) 0.26 0.23 [0.01, 0.96] NA NA NA NA NA 1.000 13653 16584
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="01_FinalModelsPlan_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw 

hypothesis(is_reactance, "exp(pressure_self_cw) > exp(pushing_self_cw)")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (exp(pressure_sel... > 0     0.81      0.95     -0.3     2.44       6.09
##   Post.Prob Star
## 1      0.86     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  #reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI'),
  
  model_rows_random = model_rows_random,
  model_rows_fixed = model_rows_fixed,
  model_rownames_random = model_rownames_random,
  model_rownames_fixed = model_rownames_fixed
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = rows_to_pack) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 4,  
      "Device-Based MVPA Log (Gaussian)" = 2, 
      "Mood Gaussian" = 2,
      #"Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2
  )
)


export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack,
  file.path("Output", paste0("AllModels", suffix, ".xlsx")), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
print(summary_all_models)
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Dichotome
exp(Est.)_hu pa_sub 95% CI_hu pa_sub exp(Est.)_nonzero pa_sub 95% CI_nonzero pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log Est. mood_gauss 95% CI mood_gauss OR is_reactance 95% CI is_reactance
Intercept 0.26*** [ 0.19, 0.35] 37.92*** [33.02, 43.54] 111.17*** [99.95, 123.93] 3.67*** [ 3.46, 3.88] 0.34*** [0.17, 0.64]
Within-Person Effects
Daily persuasion experienced 1.59*** [ 1.39, 1.85] 1.03 [ 0.98, 1.09] 1.03 [ 1.00, 1.06] 0.00 [-0.04, 0.04] 0.84 [0.68, 1.01]
Daily persuasion utilized (partner’s view) 1.34*** [ 1.19, 1.54] 1.03 [ 0.99, 1.08] 1.02 [ 0.99, 1.05] 0.02 [-0.02, 0.07] 1.12 [0.84, 1.53]
Daily pressure experienced 0.95 [ 0.69, 1.29] 0.91* [ 0.82, 1.00] 0.95 [ 0.88, 1.01] -0.03 [-0.14, 0.07] 1.97* [1.01, 4.43]
Daily pressure utilized (partner’s view) 1.49* [ 1.05, 2.29] 0.95 [ 0.86, 1.03] 0.98 [ 0.92, 1.05] -0.03 [-0.15, 0.08] 1.42 [0.58, 4.18]
Daily pushing experienced 0.94 [ 0.71, 1.28] 0.99 [ 0.92, 1.05] 1.02 [ 0.96, 1.07] 0.01 [-0.06, 0.07] 1.33* [1.03, 1.75]
Daily pushing utilized (partner’s view) 1.28* [ 1.04, 1.64] 0.96 [ 0.90, 1.02] 1.00 [ 0.96, 1.04] 0.07* [ 0.00, 0.14] 0.92 [0.62, 1.36]
Day 0.86 [ 0.64, 1.15] 0.99 [ 0.88, 1.11] 0.97 [ 0.91, 1.04] 0.26*** [ 0.15, 0.37] 1.66 [0.77, 3.60]
Own Actionplan 9.48*** [ 7.77, 11.60] 1.33*** [ 1.21, 1.45] 1.06* [ 1.01, 1.12] 0.10** [ 0.03, 0.18] 0.87 [0.46, 1.64]
Partner Actionplan 1.17 [ 0.96, 1.42] 1.08 [ 0.99, 1.17] 1.05 [ 1.00, 1.10] -0.03 [-0.11, 0.05] 0.87 [0.49, 1.54]
Daily weartime NA NA NA NA 1.00*** [ 1.00, 1.00] NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.40 [ 0.69, 2.84] 1.03 [ 0.77, 1.40] 1.10 [ 0.83, 1.46] 0.34 [-0.20, 0.90] 1.90 [0.57, 6.60]
Mean persuasion utilized (partner’s view) 1.33 [ 0.65, 2.67] 0.99 [ 0.73, 1.34] 0.98 [ 0.73, 1.30] 0.22 [-0.32, 0.78] 1.83 [0.50, 7.20]
Mean pressure experienced 0.40* [ 0.18, 0.90] 1.18 [ 0.83, 1.67] 0.99 [ 0.74, 1.33] -0.30 [-0.85, 0.24] 18.34** [2.51, 165.74]
Mean pressure utilized (partner’s view) 0.51 [ 0.22, 1.15] 0.92 [ 0.65, 1.32] 0.97 [ 0.74, 1.29] -0.31 [-0.86, 0.22] 2.29 [0.23, 19.25]
Mean pushing experienced 1.17 [ 0.41, 3.28] 1.21 [ 0.77, 1.88] 0.95 [ 0.63, 1.43] 0.21 [-0.55, 0.97] 0.86 [0.12, 6.20]
Mean pushing utilized (partner’s view) 2.07 [ 0.71, 5.84] 1.30 [ 0.82, 2.06] 1.22 [ 0.82, 1.84] 0.37 [-0.39, 1.13] 0.08* [0.01, 0.68]
Mean weartime NA NA NA NA 1.00 [ 1.00, 1.00] NA NA NA NA
Random Effects
sd(Intercept) 0.73 [0.56, 0.99] 0.30 [0.23, 0.41] 0.29 [0.23, 0.39] 0.60 [0.48, 0.78] 1.15 [0.74, 1.76]
sd(Daily persuasion experienced) 0.21 [0.03, 0.40] 0.11 [0.07, 0.17] 0.05 [0.02, 0.08] 0.04 [0.00, 0.10] 0.22 [0.01, 0.51]
sd(Daily persuasion utilized (partner’s view)) 0.19 [0.02, 0.37] 0.08 [0.04, 0.13] 0.05 [0.03, 0.09] 0.07 [0.01, 0.13] 0.49 [0.13, 0.98]
sd(Daily pressure experienced) 0.18 [0.01, 0.68] 0.06 [0.00, 0.22] 0.04 [0.00, 0.13] 0.07 [0.00, 0.24] 1.08 [0.14, 2.43]
sd(Daily pressure utilized (partner’s view)) 0.23 [0.01, 0.89] 0.06 [0.00, 0.18] 0.03 [0.00, 0.11] 0.08 [0.00, 0.26] 0.83 [0.04, 2.75]
sd(Daily pushing experienced) 0.57 [0.30, 0.98] 0.09 [0.01, 0.17] 0.07 [0.01, 0.15] 0.05 [0.00, 0.15] 0.25 [0.02, 0.62]
sd(Daily pushing utilized (partner’s view)) 0.22 [0.02, 0.54] 0.08 [0.01, 0.15] 0.03 [0.00, 0.09] 0.07 [0.00, 0.17] 0.26 [0.01, 0.96]
Additional Parameters
sigma NA NA 0.68 [0.65, 0.70] 0.57 [0.56, 0.59] 0.96 [0.94, 0.98] NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • R (version 4.4.1; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()